The rheology of dense, polydisperse granular fluids under shear. 
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The solution of the Enskog equation for the one-body velocity distribution of a moderately dense, 
arbitrary mixture of inelastic hard spheres undergoing planar shear flow is described. A generaliza- 
tion of the Grad moment method , implemented by means of a novel generating function technique, 
is used so as to avoid any assumptions concerning the size of the shear rate. The result is illustrated 
by using it to calculate the pressure, normal stresses and shear viscosity of a model polydisperse 
granular fluid in which grain size, mass and coefBcient of restitution varies amoungst the grains. The 
results are compared to a numerical solution of the Enskog equation as well as molecular dynamics 
simulations. Most bulk properties are well described by the Enskog theory and it is shown that the 
generalized moment method is more accurate than the simple (Grad) moment method. However, 
the description of the distribution of temperatures in the mixture predicted by Enskog theory does 
not compare well to simulation, even at relatively modest densities. 



Q . PACS numbers; 45.70.-n,05.20.Dd,05.60.-k,51.10.-l-y 

Introduction 

S: 

' Granular systems under rapid flow can be modelled as a fluid of inelastic hard spheres for which a variety of theoret- 
(~| ical and simulation methods may be use to explore and understand the rich phenomenology that they exhibit [1], [2], [3]. 
O Of particular interest are sheared granular flows, in which the velocity in the direction of flow varies with position 
I ^ I .. in an orthogonal direction, due to their practical relevance, accessibility to experiment and theoretical elegance. For 
a steady rate of shearing, such systems typically reach a steady state in which viscous heating, due to the shear, 
. balances collisional cooling, due to the inelastic collisions, thus giving an example of a steady state nonequilibrium 
^ ' system. A number of papers have discussed the rheology of single component sheared granular fluids [4], [5], [6], [7] in 
which all particles are mechanically identical. However, all real fluids can be expected to contain a distribution of 
particle sizes and degrees of inelasticity. The purpose of this paper is to extend these studies to dense fluids composed 
f — of an arbitrary mixture of particle sizes and inelasticities. 

The usual model for granular fluids consists of hard spheres which lose energy when they collide. Different models 
\l for the energy loss are possible and here, attention will be restricted to the simplest case in which the energy loss 
is proportional to the contribution to the kinetic energy of the normal velocities in the rest frame. This model is 
amenable to the same theoretical tools used to study elastic hard-sphere systems. In particular, it is possible to 
construct the exact Liouville equation describing the time evolution of the N-body distribution function[8],[9],[10] 
from which the Enskog and Boltzmann approximate kinetic equations follow. The latter are closed equations for 
the one-body distribution: the Enskog equation involves only the assumption of "molecular chaos" [11], [12], while 
^ , the Boltzmann equation is its low-density limit. One of the attractions of the hard-sphere models is the existence 
of the Enskog equation which allows for the description of finite density fluids outside the domain of the validity of 
the Boltzmann equation. Since one of the purposes of the present work is to provide a foundation for the study of 
>^ transport properties in realistic systems, the Enskog equation is used as a starting point. The price paid for this 
is that the results obtained must be evaluated numerically, but as discussed below, it is always possible, and quite 
$^ trivial, to take the Boltzmann limit of the final expressions and to thereby proceed analytically in this special case. 

The specific state studied here is that of uniform shear flow in which the flow is described by a velocity field 
V {r ) = a ■ r = ayx where a is the shear rate. As discussed below, this fiow admits of a uniform state with 
spatially constant density and temperature. The shear rate, temperature and degree of inelasticity are related in 
the steady state by the requirement that the viscous heating and collisional cooling balance. Although the Enskog 
equation could be solved perturbatively in the shear rate, it is difficult to carry out the expansion to sufficiently 
high order to describe physically interesting effects such as normal stresses and shear thinning (although see ref.[5] 
where this is done for the Boltzmann equation). It is instead simpler to use a moment method without making 
any assumptions about the smallness of the shear rate as has been done for elastic hard spheres [13], [14] and for the 
Boltzmann equation for sheared binary fluids [15]. A similarly method that has been used is the so-called "Generalized 
Maxwellian" of Chou and Richman[6],[7]. In fact, as shown below, these two methods can actually be viewed as special 
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cases of a generalized moment expansion about an arbitrary Gaussian state. 

An objection to using the moment method with the Enskog equation is that the calculations are technically difScult. 
In particular, the collision integrals which occur in the study of sheared fluids can be challenging to evaluate even in 
the case of the simpler Boltzmann theory (sec, e.g. ref.[15]). One contribution of this work is to present a generating 
function technique which greatly simplifies the calculations. It is shown that all integrals of interest can be obtained by 
differentiating and taking appropriate limits of a single generating function which itself simply involves the evaluation 
of a few Gaussian integrals. With this technique, it is straightforward to evaluate all results for the most general case 
of differing particle sizes and coefficients of restitution in D-dimensions. The final results in fact turn out to be as 
simple in structure as the equivalent model of sheared single-component elastic hard spheres [13], [14]. 

A question which has aroused considerable interest is the degree to which mean-field theories, like the Enskog- 
Boltzmann kinetic theory, are applicable to granular systems[16]. This provides another reason for studying the 
particular case of USF as it is possible to simulate USF with the use of modified periodic boundaries by means of the 
Lees-Edwards simulation technique[17] so as to compare the approximate kinetic theory to numerical experiments. 
Furthermore, the Enskog equation may be solved numerically by means of the Direct Simulation Monte Carlo (DSMC) 
method[18]. Both methods are used here in order to elucidate the accuracy of the analytic calculations as a method 
of solving the Enskog equation and of the Enskog approximation itself compared to simulation. 

The organization of this paper is as follows. In Section II, the Enskog theory is reviewed and the moment method 
for solving it is presented. It is shown how the moment method can be extended to allow for an arbitrary Gaussian 
reference state and the generating function formalism is introduced. The lowest non-trivial moment solutions of 
the Enskog equation are then described and used to calculate the pressure tensor for USF which thus describes the 
pressure, normal stresses and shear viscosity of the fluid in the steady state. In Section III, the solution of the moment 
equations is compared to DSMC and MD results for a model polydisperse fluid for a range of applied shear rates. The 
generalized moment method is shown to be superior to the simple (Grad) moment method and the Enskog theory 
is shown to give a good description of many rheological properties over a wide range of densities. However, the 
description of the distribution of temperatures as a function of grain size is shown to poor raising questions as to the 
accuracy of Enskog theory. The paper concludes with a discussion of the use and applicability of the results. 



THEORY 
Enskog approximation 

Consider a system of N grains which are modeled as hard spheres. Each sphere is described by a position, 'q, 
velocity, it, and a species label r. A grain of species r has mass nir while two grains of species r and s collide 
when they are separated by a distance ars ■ This array of hard-sphere diameters may be specified arbitrarily but an 
important special case is that of additive hard sphere diameters wherein each species has a fixed diameter cr^ and 
o'rs = 5 (fr + CTg). When grains i and j collide, their relative velocity after collision, it^j = 1?^ — itj is given by 

'^ij = 'Vij - (1 + "nrj Qij {Qij ■ 'Vij) (1) 

where ars is the coeflicient of restitution for collisions between grains of species r and s. The collisions are elastic 
if a^s = 1 while a^g < 1 leads to an irreversible loss of energy in each collision. Between collisions, the grains 
stream freely so that their velocities are constant. This model is a particular case of endothermic hard-sphere 
collisions in which the energy loss is proportional to the kinetic energy along the line of collision in the rest frame 
E[- — Eij = — ^1 — ctrirj^ l^rirj {Qij " ij)'^ where the reduced mass is firs = rrirms/ (m^ + mg). Finally, it is useful 
to define the momentum exchange operator for any function of the relative velocities g (^ij) as 

bijg [IJi]) = g {~v[j) = g {~Vij - (l + ar,r,) Qij (qij ■ 'vij)) (2) 
and all other velocities are left unchanged. Its inverse is 

b:^^9 (vij) = 9 ( ^ij - ^^^-^^ftj {qij ■ 'vijU . (3) 

The statistical properties of the system are determined by one and two body distribution functions, fr {l},lt;t) 
and frir2 'vi;~q2, 'v2;t) respectively. The former gives the probability density to of finding a grain of species r 
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with the given position and velocity at time t and the latter gives the joint probability for two grains. The one-body 
distribution satisfies an exact equation (the first of the Born-Bogoliubov-Green-Kirkwood-Yvon (BBGKY) hierarchy) 

+^1 • ^=r-^ /n ("51,^1; = -X] j d'q'2d'v 2T- {12) fr^r^'y'Q I, ^r-,'q 2,^ 2\t) (4) 

where the binary collision operator is 

r_ {ij) = -5 (qij - ar.rj) 

where (x) = 1 if x > and otherwise. A similar equation relates the two-body distribution to the 
three body distribution, and so on. The Enskog approximation results from noting that the combination 
S {qi2 — (Tnra) (—^12 • $12) /nrj ("9*1) 'vi',~Q2, ^21 1) picks out the pre-coUisional part of the distribution and as- 
suming that, prior to collision and at the moment of contact, the grains are uncorrelated. The specific assumption is 
that 

5 (gi2 - CTrirs) {-^12 ■ $12) /rir2 iT^ l\ ~Q2, ~V 2\t) (6) 

- 5 (gi2 - Cr^irs) © (-^12 ' 'Z12) /n , W^l i i) /rra {~q 2,~^^ 2]t) Xnn {'Ql,~Q2;t) 

where the term Xnra {~Qi^~<l2', t) .the spatial pair distribution function (pdf), accounts for spatial correlations as exist 
even in equilibrium. If it is taken to be the local-equilibrium functional of the nonequilibrium densities then the 
approximation is completely specified and is known as the Generalized Enskog Approximation[19]. 



(5) 



Hydrodynamic fields 



The local hydrodynamic fields of partial number densities nr{^,t), local velocity 'u{'q,t) , and temperature 
T Cqjt) are defined as 

nr{'q,t) = j (Tv fr{'q,l);t) (7) 
p{'q,t)li{~q,t) = X! / d'v 'mr'^fr{'q,^;t) 

r '' 

—kBT{~q,t) = X] / dPv -mrV^fr{'q,^\t) 

r ■' 

where p {~q, t) = m^rir (~q , t) is the total mass density, D is the dimensionality of the system and ks is Boltzmann's 
constant. The exact time evolution of these fields follows from Eq.(4) which gives [10] 



—rir + V • (unr) + V • j f 

^It + It ■ V + p"^ V • p" 

at 



(8) 



d_ 
dt ' 

with the number current 



n 



Dnks 



P : Vll + V ■'q 



Dnks 



7r=J (9) 
where V i = ~Ui ~ it {~qi, t), the pressure tensor P — P + P^ with kinetic and coUisional contributions 

tp^(T*,t) = Y^mr f dlTi fr{r,iri,t)ViVi (10) 

P^{'r',t) = Priir2 + otrir2) / dxidx2 qi2~q 12 {qi2 ' ^ 12^ S {qi2 - m^i^) Q {-qi2 ' it 12) 

rir2 

xfrir2{xi,X2;t) / dx S {1^ - x~q 1 - {1 - x) ~q 2) , 
Jo 
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the heat flux vector ~q = "g*^ + + ~q^^ with 

= Y.l'^- I (}^lfr{y,^l,t)ViVi^ (11) 
= ^ ^Mrurs (1 +Q;rir-2) j dXidX2 'q 12 {qi2 ■ ^ 12f S {qu - (T^h) ^ {-Ql2 ■ 1^ 12) 

x/rir2 (a;i,a;2;i) - ^ ("g • §12 y" da; J ("r - xg* i - (1 - a;) "9*2) 
'^^^{'r^,t) = V (1 - a^,r,) ^Mrurs / rfa:;irfa;2 "9i2(§i2- w'i2)^^(g'i2-o-;i;je(-gi2-^i2) 

rir2 

X /rir2 , a;2; t) / da; ^ ( - a;"g 1 - (1 - a;) "^2) 
^0 

where the center of mass velocity V = ^ " ^ ' ^'^'^ energy sink term 

C C^, i) = \^{}-- "nr.) Mriir2 / da;ida;2 12 • ^12)^ 5 (912 - cr^ijj 6 (-§"12 • V12) (12) 

rir2 

x/rir2 (a;i,X2;i)<5(^-'gi), 

Using the approximation given in Eq.(6) gives expressions for the fluxes and the heat sink which only require knowledge 
of the one-body distribution function. 

Uniform Shear Flow 

The Enskog equation is indeterminate until some boundary condition is specified. If periodic boundary conditions 

are imposed, then it is easy to see that the Enskog equation admits of a spatially homogeneous solution which 
is, however, time-dependent due to the cooling resulting from the dissipative collisions. This is the well known 
Homogeneous Cooling State (HCS). Uniform Shear Flow (USF) is another simple nonequilibrium state supported by 
this system. In USF, the density and temperature arc spatially homogeneous while the velocity field varies linearly 
with position, viz u ( r ) = a ■ r where the shear tensor, a , will be taken to be axy in a Cartesian coordinate 
system. There is therefore a flow in the x-direction which varies linearly in the y-direction. We hypothesize that 
in this case the distribution will only depend on the peculiar velocity fr {~q, l?;t) = fr(^V =1^ — 'it ■~q\t^. The 
Enskog equation then becomes 



- — Vi - a 



dt 



fr,{Vi-t) =-Y.l d^2dlf2T_{12) fr, (Vi;t) fr, (V2;t) XrrrA~QU~Q2;t) . (13) 

In fact, the linear flow field is the only one that makes the coUisional term on the right independent of position as 
follows from the observation that it will be independent of position only if li {'q'l) — li ("9*2) is a function of "9*12- 
(The function Xrir2 ("?i) "9*2; i)' evaluated in the local equiHbrium approximation, will only depend on 'qi2 if the 
density is imiform.) In fact, the requirement that "u {~qi) — it ("^2) = U ("5*12) for some field U (^12) is only 
satisfied for U {~qi2) = ~u {~q 12) (demonstrated by take ~q2 = ). One must also have that ~u {~q x) is odd, as shown 
by reversing the arguments, and that ~u = Cqi) ,for arbitrary integers p and q, as follows by iterating 

with 'q2 = —'qi- The only continuous function satisfying these constraints is one linear in 'qi, which is to say USF. 

Another consequence of the assumption of a uniform state is that the density and temperature are spatially uniform 
and it may be verified that the pressure tensor, heat-fiux vector, number flux and heating rate are all spatially uniform 
as well. The equations for the hydrodynamic fields then become 

= (14) 



dt Dnks Dnks 
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which shows that the hydrodynamic fields will also be independent of time if the viscous heating, characterized by 
the second term on the left in the temperature equation, balances the cooling described by the source term on the 
right. It is therefore consistent to hypothesize not only a spatially uniform solution to the Enskog equation, but that 
a time-independent solution exists. In this steady state, the only two quantities having the units of time are the 
temperature and the shear rate. It will therefore be the case that the relevant dimensionless control parameter is 

^* ~ ^V^^^B^"' '^'^^^^ average diameter is (a) = J^rs^r^sfTrs and the average mass is (m) = J^r^r^r- For a 
given set of material parameters, the steady state will be unique in that the value of a* as determined from eq.(14) 
will be independent of the shear rate a applied through the Lees-Edwards boundary conditions: in other words, for a 
given value of a the temperature will relax to a value such that the same value of a* is always achieved. 

Does this mean that once the fluid is moving according to the linear velocity profile, it will continue to do so forever? 
The answer is that it depends on the boundary conditions which have so far not been specified. The steady-state 
distribution hypothesized, and the sustained USF it implies, is only possible if the distribution function is compatible 
with some set of boundary conditions. For example, if the system is bounded with rigid moving walls, then the 
final distribution will depend on the detailed dynamics of collisions of grains with the wall. Here, however, we will 
assume the imposition of Lees-Edwards boundary conditions which are periodic boundaries in the co-moving frame. A 
time-independent, spatially homogeneous distribution function is in fact compatible with these boundary conditions 
and they are also amenable to application in molecular dynamics computer simulations, as discussed in more detail 
below. 

Before turning to the construction of a spatially uniform, time-independent solution of the Enskog equation, some 

comments can be made about the generality of these results. First, the only properties of the flow state used so far 
are that the flow field is a linear function of the coordinate and that the shear tensor satisfies a ■ a ~ (needed 
to convert the spatial derivative on the left in Eq.(4) into a velocity derivative in Eq.(13)). Second, the Boltzmann 
equation results from taking the lowest order term in an expansion of the integral in Eq.(13) in terms of the hard- 
sphere diameter. This results in the replacement of it 12 —» 1^12 so that, in this approximation, the collision integral 
is independent of position for any choice of the flow field. 



MOMENT SOLUTIONS TO THE ENSKOG EQUATION 

Moment Expansion 

In order to develop an approximate solution of the Enskog equation without making any restrictive assumptions 
about the size of either the shear rate or the degree of inelasticity, the distribution function is expanded in a complete 
set of polynomials about some suitable reference state [20], [21]. Here, more generally than is usually done, the reference 
state will be taken to be an arbitrary Gaussian so that the expansion takes the form 

/,(y;{nj,r) =nrdet(jy^\-^/'e^p(-V -V) 1 ^ AJ„if,„ (V2 F^'-V^ . ^^'j (15) 

< — > 

where F is a positive-definite, real, symmetric matrix and an abbreviated notation is used whereby /„ = ii...z„. 
As such, it can be written, via Cholesky decomposition, as F = F ""^/^ • ^ F ''^/^^ for some matrix F ^''-f^. The 
functions used in the expansion are the Hermite polynomials [21] given by 

Hi. (^) = (-1)" e«V2^...^e-cV2 (16) 
dci, dci„ 

so that, e.g., Hij ("?) = CiCj — Sij. They are orthonormal in D-dimensions because 

I \D/2 ^ 

_ j e-'^'I^Hi^ (^) Hj^ (^) = 5mn • (17) 

so that the coefficients of the expansion are related to the velocity moments via 

J dV fr (V; {n J , t) hi (Cr) = nrA\^ (18) 
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where Cr = T ''^/^ • V . Evaluating equation (7) relates the lower order coefficients to the hydro dynamic fields as 
nr{'q,t) = nrA" (19) 
p{'q,t)ui{'q,t) = YmrnrA''ui + YYl^rnr(y2r''^^'^y, 



implying = 1 and J2r Si f^r^r (V^ T ''^/^^ A^, = 0. The kinetic contribution to the stress tensor is 



^ 1/2 



Each species has a temperature given by 

DrirkBTr = mr J dV fr (V; {us} , y2 = 1 (mrUrTr (^(y^^) + mrUrTr (^(y^"') ^''^ ■ ■ ( r'''^) 

(21) 

Calculation of the velocity moments of Eq.(15) shows that there are actually redundant degrees of freedom in the 
sense that any change in F can always be compensated by changes in the coefficients of the expansion. These 
redundant degrees of freedom can most conveniently be eliminated by restricting the form of either F or A!j^ 
although one could imagine applying the restrictions to higher moments. There are two cases of particular interest. 
In the first, F is left unspecified and we set ^[^jj = so that all information about the second moments comes from 

the Gaussian. This will be referred to below as the Generalized Moment Expansion or GME. In the second case, F 
is specialized to a diagonal matrix by setting 

F ^ = -P— 1 (22) 
2kBTr ^ ' 

where the parameters Tr are to be determined. In this case, only the trace of the second order coefficients is set to 
zero so that 

i 

This represents an expansion about a local equilibrium state in which the temperatures of the different species are 

allowed to vary and will be referred to below as the Simple Moment Expansion or SME. Further trade-offs between 
the degrees of freedom in the reference state and those in the moment expansion are possible, but do not appear to 
offer any qualitative advantages. For example, one could use the restrictions 

- ^ (24) 



2fcsT 

r i 

SO that the reference state is simple equilibrium. This possibility will not be considered here, although there is nothing 

to rule it out in principle. Note that in no case can we force the temperatures of the subspecies to be equal since 
that would eliminate certain degrees of freedom altogether and this would lead to inconsistencies when we use the 
expansion to solve the Enskog equation. 

Substituting the general expansion given in Eq.(15) into the Enskog equation, multiplying by 

and integrating over velocities yields an infinite hierarchy of coupled equations for the moments which arc given 
explicitly in the appendix . The n-th order approximation is usually taken to consist of truncating the expansion 
in Eq.(15) to the first n terms and using the first n equations of this hierarchy to determine the moments. The 
remainder of this paper will focus on the simplest nontrivial approximations, which arc in both cases the second order 
approximation. For the GME, some simplification allows the moment approximation to be written as 

+ a (r;i-'<5nx + r--^<5,,.) = E nr^^^ (25) 
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with 



El^l^' = -27r-^ det ( "t" "t") J dVidV2dq2 Vu,Vu,T_{12) exp (-F i • 'f'^i • - • 't''^ • ^^2) (26) 
and where Xrir2 = Xr-ir2 {'Qi^'Qi + <^rsQi2',t) is independent of ~qi and gi2 in a homogeneous system. This system 

i > 

of equations suffices to determine the matrices T '^^ . For a given applied shear rate a, the steady-state temperature, 
and hence the dimcnsionlcss shear rate a*, is then determined from the scc;ond moments from Eq.(19). Exphcitly, the 
temperature is T = XrT^ where the partial temperatures are given by 



(27) 



This lowest order GME corresponds to the "Generalized Maxwellian" approximation of Chou et al[6],[7]. For the 
SME, the moment equations are 



d dlnTr 



Tr 



with 



ili2 



f~irir2 



Tjrir2 
ili2,3lh 



dt 

^ ^ ^r2Xrir2 

1"2 





dt 



(28) 



r,rir2 I 1 (r^rxr2 at^ , 1)^1^2 A^2 

lll2 ^ 9 \ tii2,3ih hh ^ »i*2,jii2 jii 



J ciy^dV 2(1-^2 H,,,, pi) r_(12)$,, (Vi;Tr,) (V2;Tr, 
J dVidV2dq2 Hi,i, pi) T_(12)$,, (y,;Tr,) Hj,j, (Ci) (VslT,,) 
dVidV2dli2 H,,i, (Ci) T_(12)$,, {Vi;Tr,) p2;T,,) H,,,^ {c 



(29) 



V2fc^7^ 



1/2 



D/2 



exp - 



2fcBT, 



The first of Eqs.(28) are a set of linear equations for the coefficients -^ili^ whereas the second can be thought of as a 
set of constraints that serve to determine the partial temperatures. 

Finally, one problem with the moment expansion in general is that the truncated distributions are not necessarily 
positive definite. An ad hoc procedure to rectiiy this problem is to resum the truncated moment expansion so that 
one writes, in the general case, 



fr (F; {n,} ,t) = nr det ( F^) tt'^/^ exp {'V ■'v' -v) ^1 + Y,^Y. ^L^/. {V2't''/^ ■ v) j 

\ n I„ J 



(30) 



where the new coefficients Aj^ are chosen so that the two series agree, term by term, up to the desired order. In 
general, the normalization constant Z must also be determined. For the second order GME, this is not an issue since 
the approximate distribution is Gaussian. For the SME, one has that 



D/2 



TT-^/^exp 



■^') ^l + ^E^U^^M2(^/2^^'■^/^.y)] (31) 



1/2 



rir det ( 1 + ^ 1 TT 



-D/2 



exp 



£7.(Y + :4')-.v') 

i\i2 / 
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which is structurally the same as the GME except that the matrix of second moments is determined through the 
linearized equations (28)-(29). Thus, in this formulation, the second order GME and SME are virtually identical 
except for the approximations used to determine the second moments. 

Generating Function for Collision Integrals 

All of the collision integrals that will be needed can be obtained from the generating function 

Z?r = -T^~'' j dVidV2d-q2 ^f](ti2),^.j expp • V^i)T_(12)exp(-yi-'r^ (32) 

by differentiating with respect to the matrices r[j and T^j and the vector A, and taking appropriate limits (such as 
A ^ 0). For example, by inspection, one has that 



Drir2 
iii2 



/< > < !■ \ 



1/2 g2 
lim 



(33) 



lim EPP 



iii2,jij2 



2kBT, 



lim 



ri r = 



^ det ( r 



1 / rrir 



lim 



lim 



-EP 



n»2,ju2 2kBTr, ksTr^ . ^ ar? 



while comparison of Eqs.(lO), evaluated for a uniform system, and (67) the coUisional contribution to the pressure is 
found to be 



1/2 



d 



m lim — Z^^^^ 

A ^ C'-'*-i2 



(34) 



The generating function is shown in appendix to be 



= ^det ( r^-)"'^'<-^ / dq I n^--% 



Z^i'-^ ( (1 + ar,r2) ^ ) - (0) 



(35) 



with 



- -a; ( A • §■) 



(36) 



X exp A • G ''^'^^ (x) • A - xw^ra A • q 
2 /"^ 

■pTi(a;) = — -1= / (w + a;)" exp (— u^) 



^rir2 = \ q- r'-i-i+ r 
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(a;) = 'f''i-i+2('f''i +'f''=) ^-2a;r'''i-i-$5-+a;%2 



rir2 — "rir2 
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and, in particular 



Fq{x) = erf (a;) - 1 

Fi{x) = ^e"^' +a;(erf(a;) - 1). 



(37) 



(The Appendix also discusses the general case of an arbitrary flow state.) The elements needed for the second order 
moment equations are worked out in Appendix where it is shown that 



'^rir2 (■'■ "I" <^''i'"2) 



triri J \Xr1r2J ^ 



(38) 



— (7 

2 '"l'"2 



Mrir2 

rrir. 



+ Fi 




from which immediately follows the coefficients for the simple moment approximation 



Drir2 
iii2 



= -2<-/(l+a,,.J 



— (7 



j dqqi. 



(1 + ar^riY 



Yr 



rir2 



Fo 



f-irir2 

jl«2 ,jlj2 



jjrir2 
n»2 ,3132 



n-1 



Til. D''! ''2 — 2rr^ 
T m »i»2,jij2 ''"rir2 
J-T2 '"^ri 



dqqi^qi^Fi 



3 irir^ 



Fi 

A*rir2 



Yr 



r\r2 



Yr, 



(1 + a^irj ^T^^rirs / dq {6i,j,qi^qj^ + Si^j,qi,qj^) Fi ( ) 

"'■ri J \ ^rir2 J 



-4(7, 



(1 + a^irs) 



rrir 



+3a,^-i(l + a,,,J^ 



/irir2 



i=0 



Yr^r2 ^ dq q^ Qi2 Qji %2 Fi 



Y 



where 



, kBTr2 

rriro 



+ 2 



(39) 



(40) 



which, together with Eqs.(28), Eq.(25), Eq.(21) and the requirement that ^^.rirTr = riT complete the specification 
of the second-order moment approximations. The coUisional contribution to the pressure is 



(41) 



The evaluation of the SME model in the Boltzmann limit, obtained by expanding in the hard-sphere diameters 
and keeping only the leading order (which here amounts to taking the limit Wrir2 and using Fq (0) = — 1 and 
-F'i(O) = ■^) is performed in Appendix . For a one component system, the GME results are in agreement with Chou 
and Richman[6],[7] while the Boltzmann limit of the SME is in agreement, for a one-component system, with the 
expressions given by Garzo[15]. In this simple case, the angular integrals can be performed analytically (see appendix 
). It is remarkable that the structure of the these terms is virtually identical to the equivalent quantities which occur 
in the elementary case of the moment solution of the Enskog equation for USE of clastic hard spheres [13], [14]. The 
practical result is that it is technically no more difficult to work with an arbitrary mixture than with a single species. 



Polydisperse Granuleir Fluids 



As an extreme example, these results can be generalized to describe a polydisperse granular fluid in which there is 
a continuous distribution of grain sizes, masses and coefficients of restitution. This is equivalent to having an infinite 
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number of species and in general one must supply the distribution of grains amongst the species, i.e. Xr as well as the 
hard-sphere diameters, a^rj and coefficients of restitution anrj- In fact, formally, the species label can be replaced 
by a continuous index over some interval, say [0, 1], and sums over species replaced by integrals over this index. In 
the event that each species has a unique hard sphere diameter (t^, and the hard-sphere diameters are additive, i.e. 



(42) 



it makes sense to replace the integrals over species labels by integrals over the distribution of hard sphere diameters. 
Specifically, the measures Xrdr become Xr^da = x {a) da where x (cr) is the fraction of grains having diameter a. 
The moment equations then become 



Ot V / nz2 \ V ^ yi2 ^ yii J 

■/ 



(43) 



n \ d(T2 X (0-2) X 



(Ti + 0-2 



(1 + a (0-1,0-2)) 



M(o"1.0"2) 

m(oi) 



Ei^i2 (^1,02) 



with 

Ei^i2 (171,02) 



(0-1,02)^1 



w (01,02) 
X (01,0-2) 

M(cri>'72) 



( r"-! (ai) • g) + ( "r-i (01) • g)^^ qi. 



(44) 



+ -(l + a(oi,a2)) ^^^^^ 



X j dqqi^qi^X^ {ai^a^) 



w (01,02) 
X (01,0-2) 



Fo 



w(0l,02) 

X (01,0-2) 



w (0-1, 02) 
X (0-1,0-2) 




w (01,02) 
X (0-1,02) 



+ 2 



X (0-1,02) 
w (0-1, 02) 



5 - (r*-i (aO + f* -1 (<T2)) 
01 + 0-2 ' 



g- a -q 



and the contributions to the pressure becomes 

i^f,, = n / doia;(oi)('r (01));' 

^n^y doido2 a: (oi) .X (02) 
X dqqi^g^X"^ (01,0-2) ^^o 



(45) 



Ol + 02 



w (01,02) 

-'^ (0-1,02) 



X I '^^ n ) (1 + a (01, 0-2)) /X (01, 02) 



_|_ 2 j ^^(0"l>'72) \ (C^l, '72) 



(0-1,0-2) 



-'^ (01,0-2) 



where some model for a (01,02) and the masses m (oi) must be supplied. The generalization of the SME expressions 
is similarly straightforward. These expressions appear formidable to implement, but if the a— integrals are performed 
using n-point Gaussian quadratures, then / dax (o) (o) J2^=i '^i^ (^j) ^ i'^i)-: where Wi are the Gaussian weights 
and the o-j are determined by the Gaussian abscissas. In this form, the calculation is identical to that for n-species 
with Xn = Wix{ai). Thus, numerically, there is no practical difference between the polydisperse fluid and a mixture 
with n-species. 



COMPARISON TO SIMULATION 
Simulation and numerical methods 



In order to evaluate the models presented here, three types of calculations were performed for three dimensional 
systems: numerical solution of the second order moment expansions, numerical solution of the Enskog equation by 
means of direct simulation monte carlo (DSMC) and molecular dynamics (MD) simulations. Comparison of the first 
two elucidates the accuracy of the second order moment approximations while comparison of both to the MD indicates 
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the accuracy of the underlying assumptions: that the state obtained is indeed USF and the assumption of molecular 
chaos. 

The focus here will be on the steady-state properties of the systems, so that when evaluating the SME and GME 

models, the time-derivatives are set identically to zero. Implementation of the (static) SME requires the numerical 
evaluation of the coefficients given in Eq.(39) and the solution of equations (28) for the partial temperatures and 
the shear rate as a function of the global temperature and coefficient of restitution (the second order moments are 
determined from Eq.(28) which are linear so that moments may be taken as given fimctions of the other parameters). 
The GME requires a similar numerical evaluation of the function E^j^'^and solution of the nonlinear moment equations, 
Eqs.(25) and (27). All numerical calculations were performed using the Gnu Scientific Library [22]. In all cases, the (two 
dimensional) numerical integrals were calculated using the GSL routine "qags" (Gauss-Kronrod 21-point integration 
rule applied adaptively until the desired accuracy is achieved) with a specification of relative accuracy of lO^^and 
absolute accuracy of 10~^for the inner integral and 10~^and 10"^ for the outer integrals. The linear equations for 
the SME moments were solved by LU decomposition and the nonlinear equations for the partial temperatures and 
shear rate were solved with the GSL "hybrids" algorithm (Powell's Hybrid method with numerical evaluation of the 
Jacobian). Convergence was considered to be achieved when the sum of the absolute value of the residuals was less 
than 10"'^. The same methods and tolerances were used to solve the nonlinear moment equations for the GME. 

The second set of calculations performed consisted of the numerical solution of the Enskog equation by means of 
the Direct Simulation Monte Carlo (DSMC) method[18]. These calculations were performed using a cubic cell with 
sides equal to the maximum of the hard sphere diameters, with 10^ points and a time step of At = O.OllTrm^t where 
Tmft is the mean free time. All calculations began from an initial configuration corresponding to the equilibrium hard 
sphere fluid. Shear flow was imposed by means of Lees-Edwards boundary conditions which are periodic boundaries in 
the Lagrangian frame[17]. For each combination of temperature, shear rates and coefficients of restitution, the initial 
configuration was relaxed over a period of lOOTmft and steady-state statistics, reported below, were then obtained by 
averaging over another lOOrmft- 

Finally, these calculations are compared below to molecular dynamics simulations. In all cases, the systems consisted 
of 500 grains and the starting configuration was the equilibrium fiuid. Shear flow was again imposed by means of 
Lees-Edwards boundary conditions. After turning on the shear flow and coUisional dissipation, the systems were 
allowed to relax for a period of 5 x 10^ collisions after which statistics were gathered for another 5 x 10^ collisions. 
Errors were computed by estimating the desired statistics using data from each period of 10^ collisions and calculating 
the standard error between the estimates (the same method was used in the DSMC calculation). In the flgures shown 
below, error bars are in general not given because in most cases, the estimated errors are comparable to or smaller than 
the size of the symbols. Exceptions to this (in the case of the temperature distributions) are explicitly commented 
upon in the text. Larger systems were not simulated as they are subject to various hydrodynamic instabilities which 
violate the assumption that the state is USF[23],[24]. 

Binary mixtures 

One cheek on the expressions given here is to compare to the results of Montanero and Garzo who have evaluated the 
SME in the Boltzmann limit for a binary mixture and compared to DSMC simulations for a variety of combinations 
of mass, diameter and density ratios. The expressions for the SME when evaluated for n (a)^ = 0.001, so as to achieve 
the Boltzmann limit, do indeed agree well with the data given in ref.[25]. A particular case is for an = o"i2 = 0-22 = 1, 
nil = 10m2, CKii = ai2 = 0^22 = 0.75 and Xi = X2 = 0.5 for which Montanero and Garzo report = —0.498 
and P^y = 0.723 from DSMC simulations and P^y = -0.498 and P^y = 0.743 from their evaluation of the SME in 
the Boltzmann limit. I flnd that the (very low density) SME gives Pj^ = -0.4981 and P^ = 0.7435 in excellent 
agreement. By comparison, the GME gives PJ^ = —0.496 and P^ = 0.726 and so is in even better agreement with 
the DSMC numerical solution of the Boltzmann equation. In addition, in the same limit, the GME is able to account 
for at least some of the normal stress differences, P^ — P^, that the SME misses (in the SME in the Boltzmann 
limit, P^ = P^) but which are clearly nonzero in the DSMC calculations. 
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Polydisperse model 

A simple model was used in which the diameters are additive, the masses scale with the diameters in the usual way 

m{a) = —pa'J^ (46) 

and the coefRcients of restitution are also additive 

a (cri,cr2) = ^ (a (cTi) + a ((72)) . (47) 
The distribution of diameters was taken to be a simple triangular distribution 




so that the average diameter is 1 and the polydispcrsity, defined as the variance divided by the average of the sizes, 
is 5 = -^u. The coefficients of restitution were assumed to scale Hnearly with the diameter with the smallest grains 
being hard, a{a = 1 — u) = 1, and the largest being soft, a (cr = 1 + m) = ao < 1, where ao is a free parameter, 
so that the average value is (a) = . The equilibrium pair structure function was evaluated using the 

approximation of ref ([26]) and the accuracy of this approximation was verified in the equilibrium (ao = 1) simulations. 
In all of the calculations reported here, the integrals over the distribution of hard-sphere diameters were performed 
using a Gauss-Lcgcndre integration scheme with 10 points. Using 5 points, the results differed by about 10%. When 
evaluating the equations for an equilibrium (ao = 1) mixture, the difference between the 5 and 10 point schemes was 
also about 10% and the absolute accuracy of the 10 point scheme compared to the known exact results was 1%. The 
MD and DSMC simulations were performed with a system obtained by a random sampling over the distribution of 
hard-sphere diameters. A variety of simulations were also performed with other samplings and it was confirmed that 
the results reported below do not vary significantly from sample to sample. 

In the following, comparisons will be made for systems at three densities: a low density fluid, n* = n {(t'^) = 0.1, a 
moderately dense fluid n* = 0.25 and a dense fluid n* = 0.5. In all cases, a value of m = 0.5 or polydispcrsity of 6 = 
20. 4% was used. All results are from a single random sampling of this distribution. For the DSMC calculations, the 
large number of points used means that the distribution is well sampled. For the MD, however, the relatively small 
number of atoms might mean that the results reported below are influenced by the particular realization used. To 
control against this, I have checked a number of data points using multiple indpendent samplings from the distribution 
and confirmed that the variation induced by different samplings is negligable, at least for the properties discussed 
below. 

Accuracy of the second moments 

Figure 1 shows the kinetic part of the stress tensor, or equivalently the second moments, as obtained from the SME, 
the GME and the DSMC. Comparison with the numerical solution of the Enskog equation, i.e. the DSMC results, 
shows that the GME gives a virtually exact estimate of the second moments at all densities and degrees of inelasticity. 
It is interesting to note that the difference between the yy and zz moments, which is zero in the Boltzmann limit 
(see Appendix ), is never very great and actually changes sign at high density. The SME is in close agreement with 
the GME. The only significant difference is in the yy and zz moments where the SME tends to underestimate the 
difference between them. Figure 2 compares the GME calculation to the MD results for the same systems. The 
calculations are in excellent agreement with the simulations at low density and remain reasonable even at the highest 
density. In particular, the xy moments are in good agreement at all densities. These results show that the GME gives 
an accurate estimate of the second velocity moments as determined by the Enskog equation and that the Enskog 
equation gives a reasonable approximation to the second moments at all densities investigated. 

Accuracy of the second moment approximation 



The next question is whether stopping at second moments is sufficient to accurately approximate the full solution 
of the Enskog equation. One measure of this is the calculation of the collisional contribution to the stress tensor. 
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FIG. 1: The absolute values of the kinetic part of the stress tensor (i.e., the second moments) normalized to nksT for three 
densities as determined by the GME (solid lines), SME (open symbols) and DSMC (filled symbols). Note that the xy moments 
are actually negative. 
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FIG. 2: Same as 1, but showing results from the GME (solid lines) and MD simulations (symbols). 



Figure 3 shows the diagonal components of this quantity as calculated from the GME and DSMC and measured in 
the MD simulations. At low density, the agreement between the GME and DSMC is good, although not quite as good 
as for the moments themselves. This shows that although higher order moments will give some small contribution, 
the GME appears, in this case at least, to be a good approximation to the solution of the Enskog equation. However, 
comparison to the MD shows the shortcomings of the Enskog equation itself. At low density, agreement is good but 
even at moderate density, considerable differences between MD and the Enskog approximation are apparent although 
the latter remains a reasonable semi-quantitative approximation. At the highest density, the differences become 
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FIG. 3: The diagonal components of the coUisional contribution to the stress tensor as a function of ao normahzed to their 
equilibrium (ao ~ 1) values. The lines are GME, the filled symbols DSMC and the open symbols MD. 



qualitative in nature. In the MD, the xx component changes non-monotonically with ao whereas the Enskog theory 
predicts a monotonic increase with increasing inelasticity. Enskog predicts little change in the yy component whereas 
in fact it drops rapidly. Only the zz component is represented at all reasonably. 



Viscoelastic properties 

Figure 4 shows the dimensionless shear rate a* as a function of ao according to the DSMC, GME and MD. All of 
these are in good agreement at all densities and values of inelasticity. This agreement is also fortunate since it means 
that any differences between Enskog and MD are not attributable to a misestimated shear rate. 

Figure 5 shows the pressure (trace of the stress tensor). In contrast to elastic hard spheres, for which the pressure 
increases with increasing shear rate[14], the pressure is nearly constant. The calculations are again all in reasonably 
good agreement with the MD. Figure 6 shows the dimensionless shear viscosity 

(49) 

a V kBT{m) ' ' 




and the viscomctric functions 

ri = i^^—iz^ (50) 



Pxx Pyy (f ) 



Pyy - Pzz {o) 



2 



which measure the normal stresses. The Enskog theory gives a very reasonable estimate for all of the viscoelastic 
properties. Although ipl is systematically underestimated, ■!/'2 the shear viscosity are well approximated at all 
densities. In all cases, the errors grow with density and decreasing qq- 



Temperature distribution 

So far, the comparisons have shown the Enskog theory and the MD to be in good agreement for bulk properties 
up to n* < 0.25. Even above this density, the physically interesting quantities - the pressure, shear viscosity and 
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FIG. 4: The reduced shear rate a* as a function of qq as determined from the GME (hues) DSMC (filled symbols) and MD 
(open symbols). 
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FIG. 5: Same as 5 but showing the reduced pressure p* = p/nksT. 



viscomctric functions, are well approximated. This picture changes when attention focuses on variations of properties 
with grain species. Figure 7 shows a comparison of the predicted temperature distribution as a function of grain 
size according to the SME, GME and DSMC for the particular value a — 0.4 as well as the zero density, Boltzmann 
limit , prediction. The SME and GME are again very good approximations to the numerical results with the former 
being slightly more accurate for the smaller grains and the latter more accurate for the larger grains for which the 
SME deviates from the Boltzmann result too slowly. The surprising result is shown in Fig. 8 which compares the 
distributions obtained from the GME and MD simulations. Although reasonable, the Enskog results are in poor 
agreement with the MD for the largest grains, especially at the lower densities and most especially for n* — 0.25. 
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FIG. 6: The reduced shear viscosity and viscometric functions, as defined in Eqs.(49)-(50) as functions of ao. The lines are 
from the GME and the symbols from MD. 

Even more surprisingly, the MD results at lower densities are in good agreement with the GME approximation to the 
Boltzmann equation. The two differences between the Boltzmann and Enskog theories are that (a) the Enskog theory 
has a higher collision frequency due to the prefactor of the pair distribution function which occurs in the collision 
term and (b) the Enskog theory accounts for the non-locality of the interactions of the grains in the collision term. It 
is hard to imagine that the second point is in error, so it seems most likely that the Enskog theory is overestimating 
the collision rate for large grains. Some support for this hypothesis comes from the fact that setting the pdf to its 
Boltzmann limit (ie. unity) increases the temperature of the largest grains by about a third of the difference between 
the Boltzmann and Enskog results for n* = 0.25. This suggests that even at low density, the Enskog theory is based 
on a poor estimate of the collision rates and so that the assumption of molecular chaos, Eq.(6), is in error. This error 
is not apparent when considering the bulk properties because the distribution of grain sizes is such that the largest 
grains make a relatively small contribution to most properties: the largest contributions come from grains near the 
middle of the distribution where the Enskog theory is relatively accurate. 

CONCLUSIONS 

In this paper, the moment approximation to the solution of the Boltzmann-Enskog kinetic theory can be generalized 
so as to represent an expansion about an arbitrary Gaussian state. This framework encompasses both the Generalized 
Maxwellian approximation as well as the simple moment expansion about local equilibrium as special cases. It shows 
in particular how corrections to the Generalized Maxwellian approximation might be calculated. 

A generating function technique was also presented as a simplified means of calculating collision integrals for the 
particular case of uniform shear flow. Although the present calculation were only performed to second order, the 
generating function technique would make higher order calculations much more feasible than more straightforward 
methods. The technique is based on the observation that the post-collisional velocities of hard spheres are linear 
functions of the pre-coUisional velocities so that pre-coUisional Gaussians remain Gaussian and integrals over such 
functions are relatively straightforward to perform. This technique is particularly valuable in anisotropic states, 
such as USE, where the usual approach to evaluating collision integrals becomes very messy. The method should be 
applicable to many other types of kinetic theory calculations. 

These general methods were applied to the particular case of arbitrary mixtures of granular fluids. It was shown, 
by comparison to DSMC simulations, that both the SME and the GME are very good approximations to the exact 
solutions to the Enskog equation for a model polydisperse granular fluid. The GME tends to be slightly more accurate 
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FIG. 7: The reduced temperature distribution T*[p) = T{a)/T as a function of grain size, a. The line is from the GME, the 
filled symbols from DSMC, the dashed line is the SME and the dotted line is from the Boltzmann equation (also in the GME 
approximation) . 
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FIG. 8: The reduced temperature distribution T'*(cr) = T{g)/T as a function of grain size, a. The full line is from the GME, 
the open symbols are from MD, and the dashed line is from the Boltzmann equation (also in the GME approximation). 



than the SME and has the additional advantage that the approximate distribution is positive definite. 

Comparison to MD simulations showed that the Enskog equation gives a good estimate of bulk properties such 
as the temperature, pressure, shear viscosity and viscomctric functions (i.e., normal stresses) over a wide range of 
coefficients of restitution and densities. Shear thinning is particularly well predicted. However, a more detailed 
examination shows that part of this agreement (particularly in the case of the viscometric functions) is due to a 
cancellation of errors while the description of the variation of temperature with grain size is in fact rather poor. The 
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fact that this agreement is so poor even at relatively low densities raises the question of whether the approximate 
kinetic theory is fundamentally lacking in some way. Possible explanations of the errors are that the local equilibrium 
pair distribution function is simply inaccurate, that the assumption of molecular chaos is violated or that the systems 
are not actually in a state of USF due, e.g., to some sort of segregation process. The exploration of these possibilities 
will be the subject of a later work. 
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MOMENT EQUATIONS 

In this appendix, the left hand side of the moment equations is developed, first for a general Gaussian state and 
then specialized to uniform shear flow. The kinetic equations take the form 

+ 1? • /. (ti, ^i; i) = Y.J [/- /^] (51) 

and the distribution is expanded as 

fr{-^,V;{ns},T)= nr det {j^) tt"^/^ exp (-W ■ ■ V^) E ^L^^- {V2T^^'^ ■ v) \ (52) 

\ra=0 ' /„ / 

where = 'v — it'^ (g). This is slightly more general than the form given in the text as we allow here for an 
arbitrary, species-dependent, linear contribution to the Gaussian. In the following, all dependence on space and time 
will not be indicated explicitly, although all quantities do in fact have such dependence. Furthermore, since we are 
only interested in the left hand side of the equation, which only involves a single species, the species label will also be 
suppressed until the end of the calculation. 

The first step is to switch variables from {~q,l?,t} to ^~q' = ~q,Ci = \/2rJj'^ {vj — uj) ,t' = f| using 
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Introducing / = det I 2 F j /, the kinetic equation becomes 



d dV 



+^ indet ( r ) + (^r J/^C. + j - indet ( r ) 
= det (2'?'-)"'^'^ J [/.,/,]. 

s 

The next step is to multiply through by Hj^ ^ and to integrate over C . These evaluations are performed using 
the basic identities, which follow directly from the definition of the Hermite polynomials, 

aH/„(c) = (c) +5/„(5,„,i?/„_, (c) (56) 



d 

da 



Hi„ {c) = SiJi„^Hi„_, {c) 



where the operator Sj^ indicates a sum over all inequivalent permutations of the indicated set of indices. Repeated 
application of these gives 

aCyHi^ (c) = Hi^^y (c) + 6,yHi^ {c) (57) 



y J y +S,i„_,Hi^_,y[c)+Syi^_,Hi^_,,(c) 

Combined with the orthonormaHty of the Hermite polynomials, and integrating by parts where needed, one then has 
that 

n-^ J dC Hi^ (c) /(c) = Ai^ (58) 

n-i j dC CyHi^ (c) -^Jic) = -S^yAi^ - SiJ,^, {Ayi^_, + 5,^_,yAi^_,) 

n-i J dC C.CyHi^ (c) ^I(c) = -SyJxy {A,i^ +SiJi^,Ai^_,) 

-Siji^x (Hi^_,.y (c)+6,yHi^_, (c)) 
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Using these, the kinetic equation becomes 

'aw'; 1 a 
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■' s 

The zeroth order equation gives 

-1 ^ -1 ^ r ^ 

so that the general equation becomes 
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Speciahzing to USF gives 



= n-idet(2'f'-) / dC-Hi^ p'^) J [/r, A] 
= .n;' J dV^Hi„ {C^) J [fr, fs] 



For the second order GME, this gives 
ar'-i/2 /^r'^i/^ 



ai 



(62) 
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Multiplying through by V^r^^ and summing over ii and 12 gives 



= 2nr-' J dVWk,Vk,^J[fr,fs]- 



For the second order SME, one has 



THE GENERATING FUNCTION 



To evaluate the various kinetic integrals, we need the generating function 

' = / dVidV2d^2 ti2. j exp (a • y 1) T_(12) exp (-Fi • T^'^i 



• V^i- V2- r'^^^- y 



0(64) 



TT-^ J dVrdV2d-q2 ^i2i j exp (- V 1 • f'''^ . y ^ - y 3 . . V2) T+{12) exp ( A • Vx) 



where the negative adjoint of the colhsion operator is 

T+ (12) = -5 {qi2 - c^nrs) (-V12 • $12) ^12 • 912 



612-1 



(65) 



and in this appendix, I continue the generalization of the first appendix and allow for an arbitrary flow state so that 
Vi = lji- 'u'^' {~qi). Using 
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= Vl-{l + ar^ri) — {Vl2- ?12)?12 
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the generating function is 



^[T = -^"''^^^2' / dVidV2dq exp (-l7i • "r-i ■Vi-V2 - r''^^ • Vl) 



(66) 



(67) 



X e {-if 12 ■ q) (V12 • g) ( exp ( A • ^ 1 - (1 + a^ra) (^12 • ?) A • gi2 ) - exp I 



(A.y,)j. 

It is enough to restrict attention to the function 

Z'-'-^x) = -2det(^'T''''T^^y^\-^ J dVidV 2 & {-If 12- q){'^ 12- q) exp (-912) 

g^2 = -Vl ■Y''' ■Vl-V2-*T''^ ■V2 + ~^ -Vl-X (^12 • ?) A • §"12. 

in terms of which the full generating function is 

= l<r.' det ( r^" r^") J dq (fl ar.r.qA [z^ ((1 + a.,.J ^) - (0)) 



(68) 



(69) 



The velocity integrals are performed by switching to relative and center of mass (CM) coordinates 



so that 



1? = V1-V2 



T? Tf Mrir2 — > 
V 2 = V V . 



In terms of the CM variables, the argument of the exponential is expanded by first using 



2 ^ 



and the remaining terms become 



.2V ■ ( ^fim'r*'-! _ ^^"f^rs ] . 



A • - X (If 12 • 912) A • q 

rriri \ ' 



where Wy^r^ = (^''^ ("?!) ~ '^^^ C^i ~ crrir25)) • (i (in USF, U'rir2 = <^rir2Q' ^ ■ The first step is to complete 
square in V 

-V ■ + 't^'A -v -2V ■( ^"''^ y " - ^fnn r'"^ . 1? + a . y 

= -(V + (Y"-' + j ■ (v + + A ■ (Y''' + Y'''^ ■ A 



with 



givmg 



512 



= -(v + Ay (j^''' + *f*''=) ■ (v + A^ + A ■ (r* + T*''^) • A 



This can be simplified by expanding the second term and using 

Mrir2 ' 



23 



so that 



512 



-(V + (Y"-' +'T'''y (v + (78) 

• r • ( r + r j • r ''^ • ^ 
+-A • (^r'^i + r'^^j -A 



Furthermore, 



V / \^ TMri mr-2 



r2 



Tir2 



so 



+ A • f r + r''M • r''^ 

+iA • ('f'''i +7'^=)"^ A 

— X ^ A ■ ^ if ■ q — XWrir2 ^ ■ 9- 



where 



B 



giving 



512 



1_> N-l ^ 

+-A • (^r''i+ r''=j -A 



(79) 



512 = -(V + Ay (Y-^ + r''^^) -(v + a) (80) 

• r '■^ • f r + r '^M • r ''^ • V 



Next, we complete the square in if using 

-1> ■ Y"-' ■ (Y"-' + Y'-'y' • f^'-i • ^ + A • (f^'-i + • f*'^^ • ^ - X (a • ^ • (81) 

= -If ■ ('f'-i-i + Y^-^-^y^ ■lf+(^- (Y"-' + • r*'-^ - a; ( A • • V 

= - + • (Y'''-^ + Y''^-'^ ■ (if + + B ■ (Y'''-^ + Y^'^-^y^ ■ B 



= -i A • (Y""^ + *f ^ • *f'''^ • (Y''^-'^ + Y'''-'^ +lx(ji ■ qj q- (Y'''-'^ + Y''^--'^ (82) 



= -(V + Ay (Y^-^ + Y"-'^ -(v + A^ (83) 

-(if+By (Y"-'-' + Y"-'-^^ ~^-(if+ 

+ B ■ (Y''^-'^ +Y''^-^^ ^ -B 



Then, using 



J exp (- (v (j''' + f'''^) -(v + A'^^dV = det (Y"-' + Y''^^ ^'"^ J exp (-F^) dV 



gives 



Z'^i''^ (a;) = -2 det f T ''^ T ''M tt^/^ det f T + T ''M 

(~B ■ ( r^'^i-i + 'r^'^-'y^ • s + i A • (7'^! + 7''=^) 



-1/2 



X exp 
X exp 



• A — XWrir2 ^ • Q 



Next, expanding 



J^^ri-i^^r2-i^ 1]^. j^^n^^rs^ ^ A = ^ A • ^'^^''^ (a;) • A 



with 



gives 



1/2 



-1/2 



u ) exp ( - A • G '"^''^ (a;) • A - ar^nra ■'^ • 9 ) • 



z*-!^^ (x) = -2 dot r ''^ r ''^ j 7r-°/2 dot r ''^ + r ''^ 

X ■ q + B ■ q — ■Wr-irg) ■ 9 ~ B ■ q + ti^nrg) 

xexp(^-^- ('f'"-i+'r''=-i) ^ 
The velocity integral is performed using 

J dfu (y—'^ ■ 9 + • q — Wrir:^ • 9 — B ■ q + Wrir2^ exp (^—^ • • 
= det (^M^ y" dli' 6 • M"^/2 -q + lB -q- w^irs) ' • M^^/^ • g - B • g + w^.^^ ) exp 



= det ( M j 7r~ 



exp (— u'^) dw' 



M-y^-q 



-1/2 _D 



M-^'^-q 



Wrir2 - B - q 



where 



so that 



2 

F„(x) = — ^ / {u' + x) exp (-w'^) du' 



Fo{x) = erf(a;)-l 
1 



Fi{x) = ^e"^ +a;(erf(a;) - 1) 
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Noting that 

-1/2 



det 



p n-l _j_ p r2 

-1/2 



1/2 



and det I ( r ''^-i + T 



-1/2 



fr'^M detfr'-=j det f r + r one has that 



The final result is then summarized as 



(92) 



with 



IdetC 



p n p 



-1/2 















V=i / 





(1 + a^irs) 



- Z'''''^ (0) 



^'■^'■^ {x) = Xr.r.Fi 



2Wr^r2 + A • r ''^-^ • gi2 _ 1 



X ( A • 



(93) 



(94) 



X exp ( - A • G ''^''^ (x) • A - a;wrir2 A • q 



G '■1''^ (x) 



^q. (f* n-i + r*r2_i j . 5- 



( (Yi) - 'u'^^ Cqi - 0-rir2?)) ' 



In this calculation, it has been implicitly assumed that F ""^ , T and , nr-2 are independent of position. However, 
this assumption is unnecessary and the same result applies for spatially dependent quantities provided the substitutions 



r-i ^ r'-i(ti) 

^-2 ^ ^-^{^^-ar.r.q) 
etc., are made and quantities involving g are brought under the integrals. 



(95) 



EVALUATION OF THE COLLISION INTEGRALS 

In this Appendix, the generating function is used to evaluate the coefficients of the moment expansions. 

Evaluation of E^L^ 



We need 



=2det(r'-i r'^^j lim 



which is evaluated using 



lim 



-Z-^-- {X) = Xr.r.F'; 



Z""^^^ ix) - Z-"^^^ (0) 



• 9 



2X, 



Xr 



(96) 



(x) (97) 



-xXr 



Wxf 



2 Mx 



Xy- 



^1^2 \ j-il I ''-'rir2 



A^rir2 / V A^ri 



(7---?)^^ ^.2 + (7---?) 



''i''2 \ 1:1' ( ^rir2 
^1 



Xf 



Xj- 



Fi 



Wr 

Xf 



q] Qii 

i2 



qiiQi 
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Substituting the explicit expression for G {x) gives 

a2 



lim 



(98) 



1 



^"j F'{ {y) + yFi' (y) + y^F^ {y) + ^-F^ {y) 



QiiQij 



so that, using F^{x) = n-F„_i(x) + (5„o2 {Fi{x) — xFo{x)), one has 



Z-^''^x,y) = -xXr,r.F^{y) 



and 



12 



--x^X,^, [yFo [y) + {2y^ + 2) F, (y)] 



(99) 



+ 2^nr-2 (l + "rir2) 



dqqi^qi^X. 



3 

rir2 



9 J 9ji 

«2 



Xr 



(100) 



Evaluation of B'-^p 



This follows by taking the appropriate limit of E^^ 
rur. 



'T-ir2 . 
i2 ' 



T>rir2 
iii2 



n r = 



lim EJT 



) 7^^'-i'-2 j dqqi,qi2Pi 



+ La^-^^il + ar,r2f''" 



rur 



dqqi^qi, 



Wr 



Y 



Fn 



Wr 



Y 

J. 



■Fi 



Wr 



Y 

± r^ 



(101) 



Y 

J r^ 



with 



Yr, 



rUr 



nir 



(102) 



Evaluation of 



We need to evaluateZ)[^^j^^^^-^^-^ = 



rriri mr 



2/sbT^i ksTr^ li™*r*. 



2fcBT: 



g prir2 
^ 1 9r'?,. n»2 
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£)rir2 

'1'2 ,Jlj2 



Then, using 



rrir, nir 



lim 



■det 



1/2 



iii2 



d 



dq Xr^r^Fl 



'4^rir2 (1+Q!rir2j 



lim 
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n^n ) ksTri ksTr. 



9ii3i2^rir2 



r2 -^^riro. <2'" rwo, 



2kBTr, 



rUr 



Fi 



Xf 



(103) 



(104) 



gives 



»1»2 Jlj2 



Using 



gives 



^^^^ I 2 '■^'■^ mr2 



lim 



Y dXr 



-Xrir2Fl 



Wr 



Xr 



--a?-l{l + ar,r2y 



lim 



j dqqi^qi2qj^qj2 



1 _i /2fcBT, 



-X^ 



Xr 



Xy 



Fi 



Xr 



Xr 



+ 2 



(I) 

6X^fl (I) 



J~)rir2 

«l«2,jlj2 



_4 D-1 /I I \ AVir2_T^-i 

rir2 V "'■ir2j ^ rir2 

X y dqqi,qi2Qji%2 (fi 



^rir2 1 


^ kBTr2 \ 


\ mr2 J 


^rir2 


\ f Wr^r2 


■^l»'2 / 


1 \^ir2 



1 l-^rir2 




/ kBTr2 \ 




1 keTr, 


V "^r2 / 



Yr^r2 j dqqi,qi2%i%2Fi(^Y^ 



rir2 
r2 



(105) 



(106) 



(107) 



Evaluation of Q^/,^.^^, 



This calculation is very similar to the preceding one. Noting that in the previous calculation we had 
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whereas from the definition 



/-■ir 11-2 
ili2,jlj2 



1 I rrir 



2 VfcsT, 



lim 



-e: 



B-^n / r * 



the present calculation will require 



rrir 



lim 



d 



we can immediately write 

(orir2 _ [ _ITT~r2_\ ( ksTr^ \ j~.rir2 

ili2,hh \kBTr2 ) \ rar, ) 



'^rir2 — '^^rir2lh1j2 



Using 



gives 



1 / rrir 



2 \kBTr, 
X / dq Fi 



\ 2 

I rr^^^ /'I J- n, "1 /fnT2_-f^ 



lim 



lim 



ri-l 



— — r 

1 Jlj2 



lim rpri r'-'^T^ 



nil 
2kBTi 



/~<rir2 
ili2,jlj2 



ksTr- 



ksTr 



rur 



-2(7^ J (1 + a^irj ^^^y^,^, / {5i^j^qi^%2 + ^i2ji9ii9j2) -P^i ( ^ 



(109) 



(110) 



(111) 



^iiii<5j2i'i (112) 



(113) 



n-l 



Evaluation of the pressure 

Recall that the collisional part of the pressure is given by 



^^i2 = I E "-i"-2Xn.2 det ( r''-^ r''-^)'^' m., Jim_ 

4 rr ^ ^ A ^ c^^i2 



Tir2 



Starting with 



gives 



1Xrir2 2 



(114) 



lim_-^Z--(x)= Xr,r2F[(^) i^^^^^-^xXr,r2Q] +Xr,r2Fli^]{-XWr,r2q^2) (US) 



lim 



d 



dAi. 



A^O '^■"■J2 

and 



= -^qi2XX^,r2 ( 'f'l' 



2 I IP' I ^'"ir2 \ . r,''^rir2 jp I '^rir2 

'x — ) ^ X — ^ I X — 



Phi2 = X]'^^i'^^2'^'^'-2^'-i^2 (1 +arir2)Mrir2 / dqQi.qi^X^^^^ (fq ( 



w. 



Ic — ] T — ^ i Ic — 



(116) 
(117) 
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SME IN THE BOLTZMANN LIMIT 



The Boltzmann limit of the coefBcients needed for the SME are 



Drir2 
iii2 



iii2)ii^2 



2- ^(l + arira) 



2 "^ri 



dqqiSi2 



-^ri ^r2 7-)rir2 



i—^r\r2 

^I.\wU2-2'^r'^r2'(l+"r,.2)^^r,r2^ / dq qi,%,^ + S,,,,%,%,J 



(118) 



-4a, 



D-l 



(1 + arir2) 



Hrir2 



rir2 







V TOr-2 V 





1- -(l + a^irj 



A*rir2 



where 



rrir 



Using the elementary integrals 



/ 



dqqiSi^ = ^kii2 



(119) 



(120) 



dqqt^qtA^qi'^ = ■^^^^{Ki2^i{i'^+ Ki'Ji{i2) 
where the area of a sphere in D dimensions is 

2^^/2 



Sd^ 



T{D/2) 



gives 



r>rir2 
iii2 



rjr\r2 
ni2,ii ^2 

11*2 i*i ^2 



(121) 

(122) 



_|_ (jrir2 (^Si-^i'^6i.^i;^ + <5i2ii ^iiij ) 
D ^ ^ (^iii2^iji2 ^iii{^i2i'2 ~^ ^iii'2^i'ii2^ 



^rir2 

C'-i'-^ = -2<-^(l+a,,„) 

£)rir2 



2 - - (1 + Q!rir-2) 



/ Hrxr2 



D-l /-I I „ ^ MnT2 



= -4<r2 (l + arir2) 



rrir 



r\r2 



1 












V TOr2 j 





2 '"ri 



5£ 



1 - - (1 +arir2) 



-1? 



£>2 +2L> 



so that 



j-)rir2 jjr2 

^1^2,^1 ^2 ^1^2 



(123) 



Qr\r2 

Zli2,l'i«2~~^'l^2 



Then, the moment equations become 



2aA^y(5ij. 



r2 



(124) 



xy 



where 



^2 

E 
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Clearly all A^- are equal fov i ^ x and the tracelessness means that A^^ = — (£>—!) Ay!y. Then 



vv 



(125) 



2aAl}y = nDB''' 
a = n D — ^ 



2 -|- 2 Ayr 



yy 



which constitute n + n + n+ l = 3n + l equations for the unknowns {^^y , ^^j/' > For a one-component 
fluid, one has that 



Y 

^ rir2 



m 



1 + 3a Sd 
1 + 3a 1 



7rm 



TTTO 



(1 + a) 



2 D + 2 

and eqs.(125) can be solved explicitly with the result that 

{l-a){D + 2) 



-2 



D V Tvm 



^yy 



3 - 3a + 2D 

Sd 



a* All = (1 - 



2^ 



(1 - a2) (1 + a) 



{D + 2) {D + l+a{D-l)) 



where 



a = a\ 



* D 

n = no . 



Recall that in this approximation 



so that 



Pyy = nksT 
Pxy = -nksT 



l + D + {D-l)a 



3 + 2£) - 3a 

1 



with 



f] = -Pxy I a = ??o 



(3 - 3a + 2D) 

1 



ID 



D + 2 



{I -a) {D + l + a{D-l)) 



3 - 3a + 2£) 



AD{D+l + a {D - 1)) 
+ 



.1-^ /-7r^ A-P + 2)r(i?/2) \ 



(126) 



(127) 



(128) 



(129) 



(130) 



(131) 
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